date.work <- "2025-04-10"
dir.create(paste0("results/", date.work), recursive = TRUE)
results.folder <- paste0("results/", date.work)
`%!in%` = Negate(`%in%`)

library(tidyverse)
library(dplyr)
library(haven)
library(labelled)
library(splines)
library(brms)
library(survey)
library(lme4)
library(ggeffects)
library(data.table)
library(glmmML)
library(srvyr)
library(writexl)
library(cmdstanr)
library(gmodels)
library(Hmisc)

options(max.print = 10000)
options(survey.adjust.domain.lonely=TRUE)
options(survey.lonely.psu="adjust")

# dynamically replace----
target.measure <- "dm_diag" #dm_diag, dm_control_diag, bp_control_diag,  statin_dm_diag, dm_control, bp_control, statin_dm, 
subset_condition <- paste0("sample_", target.measure, "==1")
AGE4 <- "FALSE"
nrow <- 100
#  Generate every 20th variable name from V1 to V2000
selected_vars <- paste0("V", seq(20, 2000, by = 20)) #V20 - V2000
## read design----
  design <- readRDS(paste0("design.", target.measure, ".WHO.rds"))

# 1. WEIGHTED DATA----
## overall----
svy_mean_overall <- function(i, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    subset_design <- subset(design, eval(parse(text=subset_condition)))

    target <- svymean(formula, na.rm = TRUE, vartype=c("var"), design=subset_design)
    return(list(target)) #  means
  })
}

results_overall_100 <- lapply(selected_vars, function(var) svy_mean_overall(var, subset_condition))
save(results_overall_100, file=paste0(results.folder, "/", target.measure, "/results_overall_100.rds"))

## country, region and income means----
svy_mean <- function(i, factor, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    target_factor <- as.formula(paste0("~", factor))
    subset_design <- subset(design, eval(parse(text=subset_condition)))
    
    target <- svyby(formula, target_factor, FUN = svymean, na.rm = TRUE, , vartype=c("var"), design=subset_design)
    return(list(target)) #  means
  })
}

## country----
results_country_100 <- lapply(selected_vars, function(var) svy_mean(var, "country_factor", subset_condition))
save(results_country_100, file=paste0(results.folder, "/", target.measure, "/results_country_100.rds"))

## region----
results_region_100 <- lapply(selected_vars, function(var) svy_mean(var, "regioncat_factor", subset_condition))
save(results_region_100, file=paste0(results.folder, "/", target.measure, "/results_region_100.rds"))

## income----
results_income_100 <- lapply(selected_vars, function(var) svy_mean(var, "countryGDPclass_factor", subset_condition))
save(results_income_100, file=paste0(results.folder, "/", target.measure, "/results_income_100.rds"))

# SEX----
svy_mean_sex <- function(i, factor, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    target_factor <- as.formula(paste0("~", factor))
    subset_design <- subset(design, eval(parse(text=subset_condition)))
    
    target <- svyby(formula, target_factor, FUN = svymean, na.rm = TRUE, design=subset_design, covmat=TRUE)
    target_d <- svycontrast(target, list(diff=c(-1,1)), add=TRUE) # male is reference
    target_r <- svycontrast(target, quote(`2`/`1`), design=subset_design) # male as reference
    
    return(list(target_d, target_r)) # difference, two means, and ratio
  })
}
results_sex_100<- lapply(selected_vars, function(var) svy_mean_sex(var, "sex_factor", subset_condition))
save(results_sex_100, file=paste0(results.folder, "/", target.measure, "/results_sex_100.rds"))

# AGE----
svy_mean_age3 <- function(i, factor, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    target_factor <- as.formula(paste0("~", factor))
    subset_design <- subset(design, eval(parse(text=subset_condition)))
    
    target <- svyby(formula, target_factor, FUN = svymean, na.rm = TRUE, design=subset_design, covmat=TRUE)
    target_d.1 <- svycontrast(target, list(diff=c(-1,1, 0)), add=TRUE) # first level as reference
    target_d.2 <- svycontrast(target, list(diff=c(-1,0, 1))) # first level as reference
    
    target_r.1 <- svycontrast(target, quote(`50-59`/`40-49`)) # first level as reference
    target_r.2 <- svycontrast(target, quote(`60-69`/`40-49`))
    return(list(target_d.1, target_d.2, target_r.1,target_r.2)) # difference, three means, and ratio
  })
}

svy_mean_age4 <- function(i, factor, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    target_factor <- as.formula(paste0("~", factor))
    subset_design <- subset(design, eval(parse(text=subset_condition)))
    
    target <- svyby(formula, target_factor, FUN = svymean, na.rm = TRUE, design=subset_design, covmat=TRUE)
    target_d.1 <- svycontrast(target, list(diff=c(-1,1, 0,0)), add=TRUE) # first level as reference
    target_d.2 <- svycontrast(target, list(diff=c(-1,0, 1,0))) # first level as reference
    target_d.3 <- svycontrast(target, list(diff=c(-1,0, 0,1))) # first level as reference
    
    target_r.1 <- svycontrast(target, quote(`40-49`/`30-39`)) # first level as reference
    target_r.2 <- svycontrast(target, quote(`50-59`/`30-39`))
    target_r.3 <- svycontrast(target, quote(`60-69`/`30-39`))
    return(list(target_d.1, target_d.2, target_d.3, target_r.1,target_r.2, target_r.3)) # difference, three means, and ratio
  })
}

if (AGE4=="TRUE") {
  results_age_100 <- lapply(selected_vars, function(var) svy_mean_age4(var, "age_cat4", subset_condition))
  save(results_age_100, file=paste0(results.folder, "/", target.measure, "/results_age_100.rds"))
} else {
  results_age_100 <- lapply(selected_vars, function(var) svy_mean_age3(var, "age_cat3", subset_condition))
  save(results_age_100, file=paste0(results.folder, "/", target.measure, "/results_age_100.rds"))
}

# BMI----
svy_mean_bmi <- function(i, factor, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    target_factor <- as.formula(paste0("~", factor))
    subset_design <- subset(design, eval(parse(text=subset_condition)))
    
    target <- svyby(formula, target_factor, FUN = svymean, na.rm = TRUE, design=subset_design, covmat=TRUE)
    target_d.1 <- svycontrast(target, list(diff=c(1,-1, 0, 0)), add=TRUE) # second level as reference
    target_d.2 <- svycontrast(target, list(diff=c(0,-1, 1, 0))) # second level as reference
    target_d.3 <- svycontrast(target, list(diff=c(0,-1, 0, 1))) # second level as reference
    target_r.1 <- svycontrast(target, quote(`<18.5`/`18.5-<25`)) # second level as reference
    target_r.2 <- svycontrast(target, quote(`25-<30`/`18.5-<25`)) 
    target_r.3 <- svycontrast(target, quote(`>=30`/`18.5-<25`)) 
    return(list(target_d.1, target_d.2, target_d.3, target_r.1,target_r.2,target_r.3)) # difference, four means and ratio
  })
}

results_bmi_100 <- lapply(selected_vars, function(var) svy_mean_bmi(var, "bmi_cat_factor", subset_condition))
save(results_bmi_100, file=paste0(results.folder, "/", target.measure, "/results_bmi_100.rds"))

# EDU----
svy_mean_edu <- function(i, factor, subset_condition) {
  suppressWarnings({
    formula <- as.formula(paste0("~", i))
    target_factor <- as.formula(paste0("~", factor))
    subset_design <- subset(design, eval(parse(text=subset_condition)))
    
    target <- svyby(formula, target_factor, FUN = svymean, na.rm = TRUE, design=subset_design, covmat=TRUE)
    target_d.1 <- svycontrast(target, list(diff=c(-1,1, 0)), add=TRUE) # first level as reference
    target_d.2 <- svycontrast(target, list(diff=c(-1,0, 1))) # first level as reference
    
    target_r.1 <- svycontrast(target, quote(`1`/`0`)) # first level as reference
    target_r.2 <- svycontrast(target, quote(`2`/`0`))
    return(list(target_d.1, target_d.2, target_r.1,target_r.2)) # difference, three means,  and ratio
  })
}

results_edu_100 <- lapply(selected_vars, function(var) svy_mean_edu(var, "educat3_factor", subset_condition))
save(results_edu_100, file=paste0(results.folder, "/", target.measure, "/results_edu_100.rds"))

# 2. EXTRACT and Summarise 100 ----
## overall----
load(paste0(results.folder, "/", target.measure, "/results_overall_100.rds"))

overall_df <- data.frame(matrix(ncol = 2, nrow = nrow))
colnames(overall_df) <- c("overall", "var_overall")

for (i in 1:nrow) {
  overall_df[i,] <- c(results_overall_100[[i]][[1]], diag(attr(results_overall_100[[i]][[1]], "var")))
}

summary_df <- overall_df %>%
  summarise(
    mean_prob = mean(overall, na.rm = TRUE),  # Average of 100 means
    within_var = mean(var_overall, na.rm = TRUE),  # Average of 100 variance
    between_var = var(overall, na.rm = TRUE)) %>%   # Variance of 100 means
  mutate(
    total_var = within_var + between_var/100,  # Total Variance
    std_err = sqrt(total_var),  # Standard Error
    lower_CI = mean_prob - 1.96 * std_err,  # Lower CI (ensuring >= 0)
    upper_CI = mean_prob + 1.96 * std_err) %>%   # Upper CI (ensuring <= 1)
  mutate(
    mean = round(mean_prob * 100, 2),
    lower = round(lower_CI * 100, 2),
    higher = round(upper_CI * 100, 2)
  ) %>%
  select(mean, lower, higher) 
writexl::write_xlsx(summary_df, paste0(results.folder, "/overall_", target.measure, ".xlsx"))

## country----
load(paste0(results.folder, "/", target.measure, "/results_country_100.rds"))
# Function to rename "var" column dynamically
rename_var_col <- function(df, index) {
  names(df)[names(df) == "var"] <- paste0("var", index)
  return(df)
}

# Apply renaming and bind
country_df <- rename_var_col(results_country_100[[1]][[1]], 1)

for (i in 2:nrow) {
  country_df <- cbind(country_df, rename_var_col(results_country_100[[i]][[1]], i))
}

# Select required columns
country_df <- country_df %>% 
  select(1, starts_with("V"), starts_with("v"))
summary_df <- country_df %>%
  mutate(
    mean = rowMeans(select(., matches("^V", ignore.case = FALSE)), na.rm = TRUE),  # Compute row means for V1 to V100
    btw_var = apply(select(., matches("^V", ignore.case = FALSE)), 1, var, na.rm = TRUE),  # Compute row variance for var1 to var100
    var = rowMeans(select(., matches("^v", ignore.case = FALSE)), na.rm = TRUE),  # Compute row means for var1 to var100
    lower = mean - 1.96 * sqrt(var + btw_var / 100),
    upper = mean + 1.96 * sqrt(var + btw_var / 100)
  ) %>% 
  mutate(mean=round(100*mean, digits=2),
         lower=round(100*lower, digits=2),
         upper=round(100*upper, digits=2)) %>% 
  select(country=country_factor, mean, lower, upper)
writexl::write_xlsx(summary_df, paste0(results.folder, "/country_", target.measure, ".xlsx"))

# region----
load(paste0(results.folder, "/", target.measure, "/results_region_100.rds"))
# Function to rename "var" column dynamically
rename_var_col <- function(df, index) {
  names(df)[names(df) == "var"] <- paste0("var", index)
  return(df)
}

region_df <- rename_var_col(results_region_100[[1]][[1]], 1)

for (i in 2:nrow) {
  region_df <- cbind(region_df, rename_var_col(results_region_100[[i]][[1]], i))
}

region_df <- region_df %>% 
  select(1, starts_with("V"), starts_with("v"))
summary_df <- region_df %>%
  mutate(
    mean = rowMeans(select(., matches("^V", ignore.case = FALSE)), na.rm = TRUE),  # Compute row means for V1 to V100
    btw_var = apply(select(., matches("^V", ignore.case = FALSE)), 1, var, na.rm = TRUE),  # Compute row variance for var1 to var100
    var = rowMeans(select(., matches("^v", ignore.case = FALSE)), na.rm = TRUE),  # Compute row means for var1 to var100
    lower = mean - 1.96 * sqrt(var + btw_var / 100),
    upper = mean + 1.96 * sqrt(var + btw_var / 100)
  ) %>% 
  mutate(mean=round(100*mean, digits=2),
         lower=round(100*lower, digits=2),
         upper=round(100*upper, digits=2)) %>% 
  select(region=regioncat_factor, mean, lower, upper)
writexl::write_xlsx(summary_df, paste0(results.folder, "/region_", target.measure, ".xlsx"))

# income----
load(paste0(results.folder, "/", target.measure, "/results_income_100.rds"))

rename_var_col <- function(df, index) {
  names(df)[names(df) == "var"] <- paste0("var", index)
  return(df)
}

# Apply renaming and bind
income_df <- rename_var_col(results_income_100[[1]][[1]], 1)

for (i in 2:nrow) {
  income_df <- cbind(income_df, rename_var_col(results_income_100[[i]][[1]], i))
}

# Select required columns
income_df <- income_df %>% 
  select(1, starts_with("V"), starts_with("v"))
summary_df <- income_df %>%
  mutate(
    mean = rowMeans(select(., matches("^V", ignore.case = FALSE)), na.rm = TRUE),  # Compute row means for V1 to V100
    btw_var = apply(select(., matches("^V", ignore.case = FALSE)), 1, var, na.rm = TRUE),  # Compute row variance for var1 to var100
    var = rowMeans(select(., matches("^v", ignore.case = FALSE)), na.rm = TRUE),  # Compute row means for var1 to var100
    lower = mean - 1.96 * sqrt(var + btw_var / 100),
    upper = mean + 1.96 * sqrt(var + btw_var / 100)
  ) %>% 
  mutate(mean=round(100*mean, digits=2),
         lower=round(100*lower, digits=2),
         upper=round(100*upper, digits=2)) %>% 
  select(income=countryGDPclass_factor, mean, lower, upper)
writexl::write_xlsx(summary_df, paste0(results.folder, "/income_", target.measure, ".xlsx"))

# SEX----
load(paste0(results.folder, "/", target.measure, "/results_sex_100.rds"))

  sex_df <- data.frame(matrix(ncol = 8, nrow = nrow))
  colnames(sex_df) <- c("diff1", "sex0", "sex1", "ratio1", 
                        "var_diff1", "var_sex0", "var_sex1",  "var_ratio1")
  
  for (i in 1:nrow) {
    sex_df[i,] <- c(results_sex_100[[i]][[1]], results_sex_100[[i]][[2]], 
                    diag(attr(results_sex_100[[i]][[1]], "var")),diag(attr(results_sex_100[[i]][[2]], "var")))
  }
  
  var.btw <- c("diff1", "sex0", "sex1", "ratio1")
  var.new <- c("diff1", "sex0", "sex1", "ratio1",
               "var_diff1", "var_sex0", "var_sex1",  "var_ratio1")
  
  summary_df <- sex_df %>%
    summarise(
      # Compute between-group variance for var.variance
      across(all_of(var.btw), \(x) var(x, na.rm = TRUE), .names = "btw_var_{.col}"),
      
      # Compute means for var.new
      across(all_of(var.new), \(x) mean(x, na.rm = TRUE))
    )
  
  # Define terms for mean estimates
  terms <- c("sex0", "sex1",  "diff1",  "ratio1")
  # Define variance columns dynamically
  var_terms <- paste0("var_", terms)      # Example: "var_sex0", "var_sex1", etc.
  btw_var_terms <- paste0("btw_var_", terms)  # Example: "btw_var_sex0", etc.
  
  # Ensure variance columns exist before computing CIs
  if (all(var_terms %in% names(summary_df)) & all(btw_var_terms %in% names(summary_df))) {
    
    # Compute final summary dataframe
    summary_df1 <- data.frame(
      term = terms,
      mean = as.numeric(summary_df %>% select(all_of(terms))),
      
      lower = as.numeric(summary_df %>% 
                           select(all_of(terms)) - 
                           1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
      ),
      
      upper = as.numeric(summary_df %>% 
                           select(all_of(terms)) + 
                           1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
      )
    )
    
  } else {
    stop("Error: Some variance columns are missing in summary_df")
  }
  
  summary_df1 <- summary_df1 %>% 
    mutate(
      mean = if_else(grepl("^ratio", term), mean, round(100 * mean, digits = 2)),
      lower = if_else(grepl("^ratio", term), lower, round(100 * lower, digits = 2)),
      upper = if_else(grepl("^ratio", term), upper, round(100 * upper, digits = 2))
    )
  writexl::write_xlsx(summary_df1, paste0(results.folder, "/sex_", target.measure, ".xlsx"))

# AGE3----
load(paste0(results.folder, "/", target.measure, "/results_age_100.rds"))
# Initialize an empty data frame
if (AGE4=="FALSE") {
age_df <- data.frame(matrix(ncol = 14, nrow = nrow))
colnames(age_df) <- c("diff1", "age0", "age1", "age2", "diff2", "ratio1", "ratio2", 
                      "var_diff1", "var_age0", "var_age1", "var_age2",  "var_diff2", "var_ratio1", "var_ratio2")

for (i in 1:nrow) {
  age_df[i,] <- c(results_age_100[[i]][[1]], results_age_100[[i]][[2]], results_age_100[[i]][[3]], results_age_100[[i]][[4]],
                  diag(attr(results_age_100[[i]][[1]], "var")),diag(attr(results_age_100[[i]][[2]], "var")), diag(attr(results_age_100[[i]][[3]], "var")),
                  diag(attr(results_age_100[[i]][[4]], "var")))
}

var.btw <- c("diff1", "age0", "age1", "age2",  "diff2",  "ratio1", "ratio2")
var.new <- c("diff1", "age0", "age1", "age2", "diff2",  "ratio1", "ratio2", 
             "var_diff1", "var_age0", "var_age1", "var_age2", "var_diff2",  "var_ratio1", "var_ratio2")

summary_df <- age_df %>%
  summarise(
    # Compute between-group variance for var.variance
    across(all_of(var.btw), \(x) var(x, na.rm = TRUE), .names = "btw_var_{.col}"),
    
    # Compute means for var.new
    across(all_of(var.new), \(x) mean(x, na.rm = TRUE))
  )

# Define terms for mean estimates
terms <- c("age0", "age1", "age2", "diff1", "diff2",  "ratio1", "ratio2")
# Define variance columns dynamically
var_terms <- paste0("var_", terms)      # Example: "var_age0", "var_age1", etc.
btw_var_terms <- paste0("btw_var_", terms)  # Example: "btw_var_age0", etc.

# Ensure variance columns exist before computing CIs
if (all(var_terms %in% names(summary_df)) & all(btw_var_terms %in% names(summary_df))) {
  
  # Compute final summary dataframe
  summary_df1 <- data.frame(
    term = terms,
    mean = as.numeric(summary_df %>% select(all_of(terms))),
    
    lower = as.numeric(summary_df %>% 
                         select(all_of(terms)) - 
                         1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
    ),
    
    upper = as.numeric(summary_df %>% 
                         select(all_of(terms)) + 
                         1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
    )
  )
  
} else {
  stop("Error: Some variance columns are missing in summary_df")
}

summary_df1 <- summary_df1 %>% 
  mutate(
    mean = if_else(grepl("^ratio", term), mean, round(100 * mean, digits = 2)),
    lower = if_else(grepl("^ratio", term), lower, round(100 * lower, digits = 2)),
    upper = if_else(grepl("^ratio", term), upper, round(100 * upper, digits = 2))
  )
writexl::write_xlsx(summary_df1, paste0(results.folder, "/age_", target.measure, ".xlsx"))
}

# AGE4----
# Initialize an empty data frame
if (AGE4=="TRUE") {
age_df <- data.frame(matrix(ncol = 20, nrow = nrow))
colnames(age_df) <- c("diff1", "age0", "age1", "age2", "age3", "diff2", "diff3", "ratio1", "ratio2", "ratio3", 
                           "var_diff1", "var_age0", "var_age1", "var_age2", "var_age3", "var_diff2", "var_diff3", "var_ratio1", "var_ratio2", "var_ratio3")

for (i in 1:nrow) {
  age_df[i,] <- c(results_age_100[[i]][[1]], results_age_100[[i]][[2]], results_age_100[[i]][[3]], results_age_100[[i]][[4]], results_age_100[[i]][[5]], results_age_100[[i]][[6]],
                       diag(attr(results_age_100[[i]][[1]], "var")),diag(attr(results_age_100[[i]][[2]], "var")), diag(attr(results_age_100[[i]][[3]], "var")),
                       diag(attr(results_age_100[[i]][[4]], "var")),diag(attr(results_age_100[[i]][[5]], "var")), diag(attr(results_age_100[[i]][[6]], "var"))
  )
}

var.btw <- c("diff1", "age0", "age1", "age2", "age3", "diff2", "diff3", "ratio1", "ratio2", "ratio3")
var.new <- c("diff1", "age0", "age1", "age2", "age3", "diff2", "diff3", "ratio1", "ratio2", "ratio3", 
             "var_diff1", "var_age0", "var_age1", "var_age2", "var_age3", "var_diff2", "var_diff3", "var_ratio1", "var_ratio2", "var_ratio3")

summary_df <- age_df %>%
  summarise( #column-wise operation
    # Compute between-group variance for var.btw
    across(all_of(var.btw), \(x) var(x, na.rm = TRUE), .names = "btw_var_{.col}"),
    
    # Compute means for var.new
    across(all_of(var.new), \(x) mean(x, na.rm = TRUE))
  )

# Define terms for mean estimates
terms <- c("age0", "age1", "age2", "age3", "diff1", "diff2", "diff3", "ratio1", "ratio2", "ratio3")
# Define variance columns dynamically
var_terms <- paste0("var_", terms)      # Example: "var_age0", "var_age1", etc.
btw_var_terms <- paste0("btw_var_", terms)  # Example: "btw_var_age0", etc.

# Ensure variance columns exist before computing CIs
if (all(var_terms %in% names(summary_df)) & all(btw_var_terms %in% names(summary_df))) {
  
  # Compute final summary dataframe
  summary_df1 <- data.frame(
    term = terms,
    mean = as.numeric(summary_df %>% select(all_of(terms))),
    
    lower = as.numeric(summary_df %>% 
                         select(all_of(terms)) - 
                         1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
    ),
    
    upper = as.numeric(summary_df %>% 
                         select(all_of(terms)) + 
                         1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
    )
  )
  
} else {
  stop("Error: Some variance columns are missing in summary_df")
}

summary_df1 <- summary_df1 %>% 
  mutate(
    mean = if_else(grepl("^ratio", term), mean, round(100 * mean, digits = 2)),
    lower = if_else(grepl("^ratio", term), lower, round(100 * lower, digits = 2)),
    upper = if_else(grepl("^ratio", term), upper, round(100 * upper, digits = 2))
  )
writexl::write_xlsx(summary_df1, paste0(results.folder, "/age_", target.measure, ".xlsx"))
}

# BMI----
load(paste0(results.folder, "/", target.measure, "/results_bmi_100.rds"))
# Initialize an empty data frame
bmi_df <- data.frame(matrix(ncol = 20, nrow = nrow))
colnames(bmi_df) <- c("diff1", "bmi0", "bmi1", "bmi2", "bmi3", "diff2", "diff3", "ratio1", "ratio2", "ratio3", 
                      "var_diff1", "var_bmi0", "var_bmi1", "var_bmi2", "var_bmi3", "var_diff2", "var_diff3", "var_ratio1", "var_ratio2", "var_ratio3")

for (i in 1:nrow) {
  bmi_df[i,] <- c(results_bmi_100[[i]][[1]], results_bmi_100[[i]][[2]], results_bmi_100[[i]][[3]], results_bmi_100[[i]][[4]], results_bmi_100[[i]][[5]], results_bmi_100[[i]][[6]],
                  diag(attr(results_bmi_100[[i]][[1]], "var")),diag(attr(results_bmi_100[[i]][[2]], "var")), diag(attr(results_bmi_100[[i]][[3]], "var")),
                  diag(attr(results_bmi_100[[i]][[4]], "var")),diag(attr(results_bmi_100[[i]][[5]], "var")), diag(attr(results_bmi_100[[i]][[6]], "var"))
  )
}

var.btw <- c("diff1", "bmi0", "bmi1", "bmi2", "bmi3", "diff2", "diff3", "ratio1", "ratio2", "ratio3")
var.new <- c("diff1", "bmi0", "bmi1", "bmi2", "bmi3", "diff2", "diff3", "ratio1", "ratio2", "ratio3", 
             "var_diff1", "var_bmi0", "var_bmi1", "var_bmi2", "var_bmi3", "var_diff2", "var_diff3", "var_ratio1", "var_ratio2", "var_ratio3")

summary_df <- bmi_df %>%
  summarise(
    # Compute between-group variance for var.variance
    across(all_of(var.btw), \(x) var(x, na.rm = TRUE), .names = "btw_var_{.col}"),
    
    # Compute means for var.new
    across(all_of(var.new), \(x) mean(x, na.rm = TRUE))
  )

# Define terms for mean estimates
terms <- c("bmi0", "bmi1", "bmi2", "bmi3", "diff1", "diff2", "diff3", "ratio1", "ratio2", "ratio3")
# Define variance columns dynamically
var_terms <- paste0("var_", terms)      # Example: "var_bmi0", "var_bmi1", etc.
btw_var_terms <- paste0("btw_var_", terms)  # Example: "btw_var_bmi0", etc.

# Ensure variance columns exist before computing CIs
if (all(var_terms %in% names(summary_df)) & all(btw_var_terms %in% names(summary_df))) {
  
  # Compute final summary dataframe
  summary_df1 <- data.frame(
    term = terms,
    mean = as.numeric(summary_df %>% select(all_of(terms))),
    
    lower = as.numeric(summary_df %>% 
                         select(all_of(terms)) - 
                         1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
    ),
    
    upper = as.numeric(summary_df %>% 
                         select(all_of(terms)) + 
                         1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
    )
  )
  
} else {
  stop("Error: Some variance columns are missing in summary_df")
}

summary_df1 <- summary_df1 %>% 
  mutate(
    mean = if_else(grepl("^ratio", term), mean, round(100 * mean, digits = 2)),
    lower = if_else(grepl("^ratio", term), lower, round(100 * lower, digits = 2)),
    upper = if_else(grepl("^ratio", term), upper, round(100 * upper, digits = 2))
  )
writexl::write_xlsx(summary_df1, paste0(results.folder, "/bmi_", target.measure, ".xlsx"))

# EDU----
load(paste0(results.folder, "/", target.measure, "/results_edu_100.rds"))
# Initialize an empty data frame
  edu_df <- data.frame(matrix(ncol = 14, nrow = nrow))
  colnames(edu_df) <- c("diff1", "edu0", "edu1", "edu2", "diff2", "ratio1", "ratio2", 
                        "var_diff1", "var_edu0", "var_edu1", "var_edu2",  "var_diff2", "var_ratio1", "var_ratio2")
  
  for (i in 1:nrow) {
    edu_df[i,] <- c(results_edu_100[[i]][[1]], results_edu_100[[i]][[2]], results_edu_100[[i]][[3]], results_edu_100[[i]][[4]],
                    diag(attr(results_edu_100[[i]][[1]], "var")),diag(attr(results_edu_100[[i]][[2]], "var")), diag(attr(results_edu_100[[i]][[3]], "var")),
                    diag(attr(results_edu_100[[i]][[4]], "var")))
  }
  
  var.btw <- c("diff1", "edu0", "edu1", "edu2",  "diff2",  "ratio1", "ratio2")
  var.new <- c("diff1", "edu0", "edu1", "edu2", "diff2",  "ratio1", "ratio2", 
               "var_diff1", "var_edu0", "var_edu1", "var_edu2", "var_diff2",  "var_ratio1", "var_ratio2")
  
  summary_df <- edu_df %>%
    summarise(
      # Compute between-group variance for var.variance
      across(all_of(var.btw), \(x) var(x, na.rm = TRUE), .names = "btw_var_{.col}"),
      
      # Compute means for var.new
      across(all_of(var.new), \(x) mean(x, na.rm = TRUE))
    )
  
  # Define terms for mean estimates
  terms <- c("edu0", "edu1", "edu2", "diff1", "diff2",  "ratio1", "ratio2")
  # Define variance columns dynamically
  var_terms <- paste0("var_", terms)      # Example: "var_edu0", "var_edu1", etc.
  btw_var_terms <- paste0("btw_var_", terms)  # Example: "btw_var_edu0", etc.
  
  # Ensure variance columns exist before computing CIs
  if (all(var_terms %in% names(summary_df)) & all(btw_var_terms %in% names(summary_df))) {
    
    # Compute final summary dataframe
    summary_df1 <- data.frame(
      term = terms,
      mean = as.numeric(summary_df %>% select(all_of(terms))),
      
      lower = as.numeric(summary_df %>% 
                           select(all_of(terms)) - 
                           1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
      ),
      
      upper = as.numeric(summary_df %>% 
                           select(all_of(terms)) + 
                           1.96 * sqrt(summary_df %>% select(all_of(var_terms)) + summary_df %>% select(all_of(btw_var_terms)) / 100)
      )
    )
    
  } else {
    stop("Error: Some variance columns are missing in summary_df")
  }
  
  summary_df1 <- summary_df1 %>% 
    mutate(
      mean = if_else(grepl("^ratio", term), mean, round(100 * mean, digits = 2)),
      lower = if_else(grepl("^ratio", term), lower, round(100 * lower, digits = 2)),
      upper = if_else(grepl("^ratio", term), upper, round(100 * upper, digits = 2))
    )
  writexl::write_xlsx(summary_df1, paste0(results.folder, "/edu_", target.measure, ".xlsx"))